Visualize integration (level 2): GCBC
2 junio 2021
1 Introduction
The objectives of this notebook are:
Visualize how well we could remove technical variability associated with the assay (scATAC or multiome). We should see a high degree of intermixing between the 2 assays.
Plot several markers of potential doublets or problematic cells: - scRNAseq doublets determined by 10X Multiome - Scrublet doublet predictions - nCount and nFeatures of scATAC - Chromatin Signature - Projection of the scRNA-seq clusters defined at level 3
The broader objective of level 2 is to eliminate most remaining doublets and poor-quality cells, as we will discuss in future notebooks. These visualizations will allow us to explore this. We will also plot their location in the UMAP obtained at level 1.
2 Pre-processing
2.1 Load packages
library(Seurat)
library(tidyverse)
library(reshape2)
library(ggpubr)2.2 Parameters
# Paths
path_to_obj <- str_c(
here::here("scATAC-seq/results/R_objects/level_2/"),
params$cell_type,
"/",
params$cell_type,
"_integrated_level_2.rds",
sep = ""
)
path_to_obj_RNA <- str_c(
here::here("scRNA-seq/3-clustering/3-level_3/"),
params$cell_type,
"_clustered_level_3_with_pre_freeze.rds",
sep = ""
)
path_to_doublets <- here::here("scRNA-seq/3-clustering/2-level_2/tmp/doublets_multiome_df_all.rds")
# Functions
source(here::here("scRNA-seq/bin/utils.R"))2.3 Load data
# Seurat object
seurat <- readRDS(path_to_obj)
seurat## An object of class Seurat
## 270784 features across 27497 samples within 1 assay
## Active assay: peaks_macs (270784 features, 267464 variable features)
## 3 dimensional reductions calculated: umap, lsi, harmonyseurat_RNA <- readRDS(path_to_obj_RNA)
seurat_RNA## An object of class Seurat
## 37378 features across 81653 samples within 1 assay
## Active assay: RNA (37378 features, 0 variable features)
## 3 dimensional reductions calculated: pca, umap, harmonytonsil_RNA_annotation <- seurat_RNA@meta.data %>%
rownames_to_column(var = "cell_barcode") %>%
dplyr::filter(assay == "multiome") %>%
dplyr::select("cell_barcode", "seurat_clusters")2.3.1 Visualization of the data
p1 <- DimPlot(seurat,
pt.size = 0.1)
p2 <- DimPlot(seurat_RNA,
group.by = "seurat_clusters",
pt.size = 0.1, label = T) + NoLegend()
p1 + p2
3 Assess integration
3.1 Technical variability |
p_assay <- plot_split_umap(seurat, var = "assay")
p_assay
3.2 Technical variability ||
p_assay <- plot_split_umap(seurat, var = "hospital")
p_assay
3.3 Biological variability |
p_assay <- plot_split_umap(seurat, var = "age_group")
p_assay
3.4 Biological variability ||
p_assay <- plot_split_umap(seurat, var = "sex")
p_assay
4 Spot potential doublets
4.1 Multiome scRNA-seq data
4.1.1 Doublets detected in scRNA-seq data from Multiome.
Here, we can see the count of doublets detected by cell-type. Note that for the level1 annotation there are not doublets detected in the Naive and Memory B cell cluster.
multiome_doublets <- readRDS(path_to_doublets)
dfm = melt(table(multiome_doublets$cell_type))
dfm$value = as.numeric(as.character(dfm$value))
ggbarplot(dfm, x = "Var1", y = "value",
fill = "Var1",
palette = "jco",
sort.val = "desc",
sort.by.groups = FALSE,
x.text.angle = 90
)
doublets_cells <- colnames(seurat)[which(colnames(seurat) %in% multiome_doublets$barcode)]
length(doublets_cells)## [1] 15714.2 Visualize the projection of the doublets cells onto the scATAC-seq UMAP
DimPlot(
seurat, reduction = "umap",
cols.highlight = "darkred", cols= "grey",
cells.highlight= doublets_cells,
pt.size = 0.1
)
## Scrublet prediction
# Scrublet
DimPlot(seurat, group.by = "scrublet_predicted_doublet_atac")
table(seurat$scrublet_predicted_doublet_atac)##
## FALSE TRUE
## 24936 25614.3 QC metrics
qc_vars <- c(
"nCount_peaks",
"nFeature_peaks"
)
qc_gg <- purrr::map(qc_vars, function(x) {
p <- FeaturePlot(seurat, features = x)
p
})
qc_gg## [[1]]
##
## [[2]]
4.4 Annotation probability
qc_vars <- c(
"annotation_prob")
qc_gg <- purrr::map(qc_vars, function(x) {
p <- FeaturePlot(seurat, features = x)
p
})
qc_gg## [[1]]
4.5 Chromatin Signature
qc_vars <- c("NBC.MBC", "GCBC", "PC", "CD4.T", "Cytotoxic",
"myeloid", "FDC", "PDC")
qc_gg <- purrr::map(qc_vars, function(x) {
p <- FeaturePlot(seurat, feature = x, max.cutoff = 4, min.cutoff = -4) + scale_color_viridis_c(option = "magma")
p
})
qc_gg## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]